function [ rf_t, rh_t, rrf_t, AP_S, mapp, s_growth_M, tau_t, s_agg_c ] = lsap_zlb_duration( filename )
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

%----------------------------------------------------------------------------
% Some Parameters
%----------------------------------------------------------------------------
T = 301;    %Number of periods 
var_foc = 29;      % Number of variables in one period
pos_value = 3;     % Position of var_value in A
truc = (-1):(T-2);         % Horizontal axis when drawing with level
spline_truc = (-1):(1/4):(T-2);     %Horizontal axis when drawing with spline 
size_struc = 1201;  % size of spline truc
deltab = 0.5;

%----------------------------------------------------------------------------
% Mapping name of variables to its position
%----------------------------------------------------------------------------
mapp = containers.Map;
mapp('gamma')=1; mapp('cb')=2; mapp('Rf')=3; mapp('Rm')=4; mapp('ql')=5;
mapp('mur')=6; mapp('muc')=7; mapp('ip')=8; mapp('pie')=9; mapp('n')=10;
mapp('m')=11; mapp('bh')=12; mapp('k')=13; mapp('tau')=14; mapp('s')=15;
mapp('ch')=16; mapp('etaz')=17; mapp('lama')=18; mapp('lamb')=19; mapp('etab')=20;
mapp('i')=21; mapp('pm')=22; mapp('y')=23; mapp('l')=24; mapp('u')=25;
mapp('Rn')=26; mapp('x')=27; mapp('up')=28; map('w')=29;

%----------------------------------------------------------------------------
%Read data file into matrix A (time X var_position X var_value) 
%----------------------------------------------------------------------------
A = zeros(T*var_foc, pos_value);
A = dlmread(filename);

%----------------------------------------------------------------------------
%Create matrix AT: row (time)X col (var_position)
%----------------------------------------------------------------------------
AT = zeros(T, var_foc);
for t=1:T
    for i=1:var_foc
        AT(t,i)= A((t-1)*var_foc+i, pos_value);
    end;
end;

%----------------------------------------------------------------------------
%Create matrix AP: Percentage deviation from steady state
%AP contains 2 columns: time X percentage_deviation
%----------------------------------------------------------------------------
% Vector stores the steady state "ss"
ss = zeros(1, var_foc);
for i=1:var_foc
    ss(1,i) = AT(T,i);
end;

% Matrix AP
AP = zeros(T, var_foc);
for t=1:T
    for i=1:var_foc
        AP(t,i)= (AT(t,i)-ss(i))/ss(i)*100;
    end;
end;

% Variable tau=0 and x_ss=0, so we cannot divide by 0, fix AP at 'tau' and 'x' position
for t=1:T
    AP(t,mapp('tau'))=0;
    AP(t,mapp('x'))=0;
end;
        


%----------------------------------------------------------------------------
%Graph 1: Federal Funds rate (level0, real FFR, real lending rate (anually)
%---------------------------------------------------------------------------- 
rf_t = zeros(1,T);   % Federal Funds Rate Level
rrf_t = zeros(1,T);  % Real Federal Funds Rate
rh_t = zeros(1,T);   % Real Borrowing rate 
for t=1:T
    rf_t(t) = (AT(t,mapp('Rf'))-1)*400;
end;
for t=1:(T-1)
    rrf_t(t) = (AT(t,mapp('Rf'))/AT(t+1,mapp('pie'))-1)*400;
    rh_t(t) = ((deltab+ deltab*AT(t+1,mapp('ql')))/AT(t,mapp('ql'))/AT(t,mapp('pie'))-1)*400;
end;
rrf_t(T) = rrf_t(T-1);
rh_t(T) = rh_t(T-1);
rrf_t(1) = rrf_t(T);   % Before the unexpected shock
rh_t(1) = rh_t(T);

%----------------------------------------------------------------------------
%Graph 2: C_bankers vs C_household vs Aggregate C
%----------------------------------------------------------------------------
agg_c = zeros(T, 1);   %Aggregate consumption
percent_agg_c = zeros(T, 1); % Percentage deviation of aggregate consumption from steady state

for t=1:T
    agg_c(t) = AT(t,mapp('cb'))+ AT(t,mapp('ch'));
end;
agg_c_ss = agg_c(T); % Steady state of aggregate consumption
for t=1:T
    percent_agg_c(t) = (agg_c(t)-agg_c(T))/agg_c(T)*100;
end;
s_agg_c = spline(truc, percent_agg_c, spline_truc);

%----------------------------------------------------------------------------
%Graph 4: Nominal growth M(t)/M(t-1), N(t)/N(t-1), B(t)/B(t-1)
%----------------------------------------------------------------------------
growth_M = zeros(1,T);
growth_B = zeros(1,T);
growth_N = zeros(1,T);
for t=2:T
    growth_M(t) = (AT(t,mapp('m'))/AT(t-1,mapp('m'))*AT(t,mapp('pie'))-1)*100;
    growth_B(t) = (AT(t,mapp('bh'))/AT(t-1,mapp('bh'))*AT(t,mapp('pie'))-1)*100;
    grwoth_N(t) = (AT(t,mapp('n'))/AT(t-1,mapp('n'))*AT(t,mapp('pie'))-1)*100;
end;

s_growth_M = spline(truc, growth_M, spline_truc);
s_growth_B = spline(truc, growth_B, spline_truc);
s_growth_N = spline(truc, growth_N, spline_truc);

%----------------------------------------------------------------------------
%Create matrix AP_S: Percentage deviation from steady state (spline)
%AP contains 2 columns: time X percentage_deviation
%----------------------------------------------------------------------------       
AP_S = zeros(var_foc, size_struc);
for i=1:var_foc
    AP_S(i,:)= spline(truc, AP(:,i), spline_truc);
end;        

%----------------------------------------------------------------------------
%Helicopter drop
%----------------------------------------------------------------------------
tau_t = zeros(1, T);
for t=1:T
    tau_t(t)= AT(t, mapp('tau'));
end;














end

